CopulaHiC is R package for differential analysis of Hi-C data. It takes advantage of significant correlations of main diagonals between different Hi-C data sets (cell lines, experiments, etc.) - usually first 200 to 300 diagonals from main diagonal are considered. CopulaHiC uses copulas to model these dependancies and then quantifies deviations from the model in a probabilistic way.

This vignette contains description of our method: - basic introduction on Hi-C data and its related biases, - issues with Hi-C differential analysis, - joint normalization of Hi-C data via copula models, - detection of differential significantly interacting regions.

If you want to spare yourself extra reading and would like to jump straight to analysis part see quick start chapter. It contains step by step introduction to Hi-C differential analysis with our package on example Hi-C data set.

1 Quick Start

Load package:

library("CopulaHiC")

To compare Hi-C experiments you need 2 files with Hi-C matrices. In this tutorial we will use Hi-C data sets of human IMR90 and human MSC cells lines from GSE63525 and GSE52457 studies respectively available as sample data in CopulaHiC package. For simplicity we will only use chromosome 18. First read data:

mtx.fname.imr90 <- system.file("extdata", "IMR90-MboI-1_40kb-raw.npz",
                               package = "CopulaHiC", mustWork = TRUE)
mtx.fname.msc <- system.file("extdata", "MSC-HindIII-1_40kb-raw.npz",
                             package = "CopulaHiC", mustWork = TRUE)

hic.comparator <- HiCcomparator(mtx.fname.imr90, mtx.fname.msc, mtx.names = c("18"), do.pca = TRUE)

1.1 Contact maps and A/B compartments

Illustrate Hi-C contact maps with A/B compartments:

dense.imr90 <- sparse2dense(hic.comparator$maps1[["18"]][c("i","j","val")],
                            N = hic.comparator$maps.dims[["18"]][1,"n.rows"])
dense.msc <- sparse2dense(hic.comparator$maps2[["18"]][c("i","j","val")],
                          N = hic.comparator$maps.dims[["18"]][2,"n.rows"])

# there are some troubles with adjusting margins, for a momement it must be done manually
layout(matrix(c(1,2,3,4), ncol = 2, byrow = TRUE), widths = c(2,2), heights = c(2,1))
par(mar = c(2,1,0.5,1))
plot_contact_map(dense.imr90)
par(mar = c(1,1,1,1))
plot_contact_map(dense.msc)
plot_pc_vector(hic.comparator$pc1.maps1[["18"]])
plot_pc_vector(hic.comparator$pc1.maps2[["18"]])

IMR90 and MSC contact maps with A/B compartments of human chromosome 18.

1.2 Contact maps and TADs

Determine TADs for both Hi-C maps and illustrate them:

tads.imr90 <- map2tads(dense.imr90, resolution = 40*1000, window.bp = 400*1000, delta.bp = 80*1000)
tads.msc <- map2tads(dense.msc, resolution = 40*1000, window.bp = 400*1000, delta.bp = 80*1000)

plot_with_inset(list(dense.imr90), c(600,900), c(600,900), args.regions = list(tads.imr90))
plot_with_inset(list(dense.msc), c(600,900), c(600,900), args.regions = list(tads.msc))

IMR90 contact map of human chromosome 18 with TADs.

MSC contact map of human chromosome 18 with TADs.

1.3 Coverage and decays

Compare coverage and decays for IMR90 and MSC, for decays use log10 scales for both X and Y axis:

library("ggplot2")

coverage <- dominating_signal(hic.comparator, which.signal = "coverage")

ggplot(coverage, aes(x = i, y = sum.contacts, color = dataset)) +
  geom_point(size = 0.5) +
  geom_smooth(alpha = 0.5) +
  facet_wrap(~ name, ncol = 1, scales = "free") +
  theme(legend.position = "bottom")

decay <- dominating_signal(hic.comparator, which.signal = "decay")

ggplot(decay[decay$diagonal != 0,], aes(x = diagonal, y = mean.contacts, color = dataset)) +
  geom_point(size = 0.5) +
  scale_x_log10() +
  scale_y_log10() +
  facet_wrap(~ name, ncol = 1, scales = "free") +
  theme(legend.position = "bottom")

Coverage comparison between IMR90 and MSC Hi-C data of human chromosome 18.

Decay comparison between IMR90 and MSC Hi-C data of human chromosome 18.

1.4 Fold Change and Differential maps

Calculate fold change and differential maps IMR90 / MSC and IMR90 - MSC respectively. Illustrate results with zoomin of 600 - 900 bins region. By default when plotting fold change maps log10 scale is used. When plotting differential map one should indicate sqrt transformation of data for better visibility (see examples below):

merged <- merge(hic.comparator)
merged.18 <- merged[["18"]]

merged.18$fc <- merged.18$val.x / merged.18$val.y
merged.18$difference <- merged.18$val.x - merged.18$val.y

dense.fc <- sparse2dense(merged.18[c("i","j","fc")], N = hic.comparator$maps.dims[["18"]][1,"n.rows"])
dense.diff <- sparse2dense(merged.18[c("i","j","difference")],
                           N = hic.comparator$maps.dims[["18"]][1,"n.rows"])

plot_with_inset(list(dense.fc, fc = TRUE, colors.pal = c("blue","white","red")), c(600,900), c(600,900), which.map = "contact.map")
plot_with_inset(list(dense.diff, breaks = 100, sqrt.transform = TRUE), c(600,900), c(600,900), which.map = "diff.map")

Fold Change map of IMR90 / MSC of human chromosome 18.

Differential map of IMR90 - MSC of human chromosome 18.

1.5 Diagonal correlations

Check correlations between corresponding diagonals of IMR90 and MSC:

library("reshape2")

decay.cors <- decay_correlation(hic.comparator)
# wide to long
decay.cors.long <- reshape2::melt(decay.cors[c("name","diagonal","pcc","rho","tau")],
                                  id.vars = c("name","diagonal"), variable.name = "correlation",
                                  value.name = "coefficient")
# remove 0 diagonal (as it is non informative anyways) and illustrate results
ggplot(decay.cors.long[decay.cors.long$diagonal != 0,],
       aes(x = diagonal, y = coefficient, color = correlation)) +
  geom_point(size = 0.7) +
  scale_x_continuous(breaks = c(1,seq(0, max(decay.cors.long$diagonal), by = 250)[-1])) +
  facet_wrap(~ name, ncol = 1, scales = "free") +
  theme(legend.position = "bottom")

dc <- decay.cors[c("name","diagonal","pearson.pval","spearman.pval","kendall.pval")]
colnames(dc) <- c("name","diagonal","pearson","spearman","kendall")
decay.sig.long <- reshape2::melt(dc, id.vars = c("name","diagonal"),
                                 variable.name = "correlation", value.name = "pval")
decay.sig.long$neg.log.pval <- -log10(decay.sig.long$pval)

ggplot(decay.sig.long[decay.sig.long$diagonal != 0,],
       aes(x = diagonal, y = neg.log.pval, color = correlation)) +
  geom_point(size = 0.7) +
  geom_hline(yintercept = -log10(0.01)) +
  annotate("text", 1800, -log10(0.01), vjust = -1, label = "pval = 0.01") +
  scale_x_continuous(breaks = c(1,seq(0, max(decay.sig.long$diagonal), by = 250)[-1])) +
  ylab("-log10(pval)") +
  facet_wrap(~ name, ncol = 1, scales = "free") +
  theme(legend.position = "bottom")

Correlations between diagonals of IMR90 and MSC Hi-C data of human chromosome 18.

Significances of correlations between diagonals of IMR90 and MSC Hi-C data of human chromosome 18.

1.6 Copula modelling

Construct Hi-C copula models (if you are running this on a cluster you can increase number of cores - it will speed up model construction significantly):

hic.copula <- HiCcopula(hic.comparator, diagonals = 1:240, n.cores = 1)

1.6.1 Marginal models and copula model

Plot marginals gamma fit for diagonal 10 and copula fit for diagonals, say 1, 10, 50 and 200:

model.18.1 <- hic.copula$model[["18"]][["1"]]
model.18.10 <- hic.copula$model[["18"]][["10"]]
model.18.50 <- hic.copula$model[["18"]][["50"]]
model.18.200 <- hic.copula$model[["18"]][["200"]]

library("fitdistrplus")
library("VineCopula")

plot(model.18.10$marginal.x)
plot(model.18.10$marginal.y)

c1 <- plot(model.18.1$bf.copula, main = "diagonal 1")
c2 <- plot(model.18.10$bf.copula, main = "diagonal 10")
c3 <- plot(model.18.50$bf.copula, main = "diagonal 50")
c4 <- plot(model.18.200$bf.copula, main = "diagonal 200")

grid.arrange(c1,c2,c3,c4, ncol=2)

Gamma fit of marginal distribution X (IMR90) of diagonal 10.

Gamma fit of marginal distribution Y (MSC) of diagonal 10.

Copula fit of U (IMR90) vs V (MSC) for diagonals 1,10,50 and 200.

1.6.2 Differential (significance) matrix

When background model is constructed differential map can be computed (and visualized):

diffmap <- hicdiff(hic.copula)
diffmap.dense <- hicdiff2mtx(diffmap, hic.copula$maps.dims)

# plot diffmap
diffmap.18.dense <- diffmap.dense[["18"]]
plot_with_inset(list(diffmap.18.dense, breaks = 100), c(1300,1700), c(1300,1700), which.map = "diff.map")

Differential/Significance map of human chromosome 18 IMR90 vs MSC cell lines.

1.6.3 Detection of Long Range differential interactions

Finally groups of cells with significant depletion/enrichment can be identified (with some precision) as rectangle like regions containing connected components:

lrdi <- differential_interactions(hic.copula)

# get rectangle like regions
regions <- lrdi[["18"]]$interacting.regions
# plot them --> only those which contain at least 5 cells inside connected component
plot_with_inset(list(diffmap.18.dense, breaks = 100), c(1300,1700), c(1300,1700), which.map = "diff.map",
                args.regions = list(regions[regions$n.cells >= 5, 2:6], pal.colors = c("blue","red")))
plot_with_inset(list(diffmap.18.dense, breaks = 100), c(900,1200), c(900,1200), which.map = "diff.map",
                args.regions = list(regions[regions$n.cells >= 5, 2:6], pal.colors = c("blue","red")))

Detection of differential long range interactions in IMR90 vs MSC map. Bilinear model fit to depleted (negative) significances in order to determine threshold for long range interactions.

Detection of differential long range interactions in IMR90 vs MSC map. Bilinear model fit to enriched (positive) significances in order to determine threshold for long range interactions.

Differential/Significance map of human chromosome 18 IMR90 vs MSC cell lines with long range differential interactions (zoomin 1).

Differential/Significance map of human chromosome 18 IMR90 vs MSC cell lines with long range differential interactions (zoomin 2).

2 About Hi-C protocol

Hi-C is the high resolution variant of 3C (Chromosome Conformation Capture) techniques, which samples chromatin interactions genomewide (Lieberman-Aiden et al. 2009). The result of Hi-C experiment is library of paired-end reads indicating pairs of chromatin regions that were in close 3D proximity in initial sample (i.e. interacting regions). This library is then mapped to reference genome to locate interacting regions and interactions are counted and summarized in the form of non negative matrices called contact maps. To construct contact map one must first specify resolution parameter - number of basepairs per unit (typically between mega- or kilobasepairs). Correct selection of resolution parameter was studied in literature and there are some guidelines as to how approach this problem (Lajoie, Dekker, and Kaplan 2015; Yang et al. 2017). A contact map cell with coordinates i,j contains number of interactions between genomic regions i and j inside the sample captured by the protocol. An important note here is that in most cases Hi-C is performed on sample containing billions of cells (rather than single, isolated cell) and therefore resulting contact maps contain sum of interactions over (possibly) distinct cell populations. There are 2 types of contact maps: inter- and intrachromosome (symmetric). Usually due to effect called chromosome territories of interest are only intrachromsome matrices.

Contact map of human chromosome 18 IMR90 cell line in 40kb resolution.

Hi-C is a complex protocol with many sources of bias (Yaffe and Tanay 2011), which make analysis incorporating this type of data challenging. Therefore almost any biological study utilizing Hi-C data performed without proper normalization will lead to nonmeaningfull results.

Here we consider a problem of comparing 2 Hi-C datasets. Depending on experimental design it may be for example different cell lines or different experimental conditions. This sort of studies apper relatively often in literature indicating high importance of this problem (Le Dily et al. 2014; Dixon et al. 2015; Lun and Smyth 2015). We present key issues arising in this type of assays and suggest a model, which have potential to improve the reliability of Hi-C differential analysis.

3 Hi-C differential analysis

Hi-C differential analysis aims to discover and quantify regions which are significantly depleted or enriched in interactions between 2 contact maps of interest. Due to large between sample variability of coverage and contact decay dominating Hi-C signal intensity, direct comparisons of contact maps would usually end up with misleading conclusions.

Coverages (sum of contacts per region) of human chromosome 19 from 2 Hi-C experiments: IMR90 and MSC cell lines.

Decays (mean of contacts per diagonal) of human chromosome 19 from 2 Hi-C experiments: IMR90 and MSC cell lines.

As an example imagine a scenario where in first Hi-C experiment sequencing depth is 3 times that of the second experiment. After conducting naive comparison we are likely to conclude that almost any pair of regions is enriched in interactions in experiment 1 with respect to 2.

Fold Change map of IMR90 / MSC of human chromosome 19.

Although there are some techniques that take decent care of coverage (Imakaev et al. 2012), strong contact decay signal can still persist and limit ability to detect local differences in contact intensity.

Coverages (sum of contacts per region) of human chromosome 19 from 2 Hi-C experiments: IMR90 and MSC cell lines after Iterative Correction.

Decays (mean of contacts per diagonal) of human chromosome 19 from 2 Hi-C experiments: IMR90 and MSC cell lines after Iterative Correction.

Fold Change map of IMR90 / MSC of human chromosome 19 with both maps iteratively corrected.

A possible workaround could be joint normalization method, which tries to model similarities in coverage and contact decay between both Hi-C experiments simultaneously and then seeks for deviations from the model.

4 Copula modelling of Hi-C data

Our method is based on an observation that the pattern of interactions between corresponding diagonals of Hi-C maps is similar and highly significant for pairs of regions located nearby as measured by linear chromosomal distance. As an example consider below images illustrating Pearson, Spearman and Kendall correlation between pairs of consecutive diagonals from pair of contact maps. To calculate correlations only cells containing interactions in both experiments (contact maps) are used.

Correlations between corresponding diagonals (i.e. 1 vs 1, 2 vs 2, etc) of IMR90 replicate 1 vs replicate 2 Hi-C contact map of human chromosome 19.

Significances of the above correlations.

Altough correlations between experiments originating from different cell lines are much lower they are still highly significant indicating of high similarities of the interactions pattern.

Correlations between corresponding diagonals (i.e. 1 vs 1, 2 vs 2, etc) of IMR90 vs MSC Hi-C contact map of human chromosome 19.

Significances of diagonals correlations.

The trends shown above are representative for many other cell types and chromosomes analyzed by us.

The above observation prompted us to model dependence between Hi-C contact maps as global “uninteresting” signal and then try to find deviations from the model as local, relevant differences. The common way to model dependencies between datasets is with the use of copulas. Copulas are multivariate distributions with uniform marginals. They allow to separate modelling of marginal distributions from modelling of dependancies between random variables. Copulas have a long history in data modelling, for more details on the topic see for example (embrechts2001modelling; Trivedi, Zimmer, and others 2007). A great advantage of using copulas is that there exist rich and well documented R packages dedicated to copula modelling, which makes model selection quick and straightforward process (Hofert et al. 2017; Jun Yan 2007; Ivan Kojadinovic and Jun Yan 2010; Marius Hofert and Martin Mächler 2011; Schepsmeier et al. 2018).

Our package is build on top of VineCopula library, which can easily handle model selection for bivariate distributions. The pipeline can be summarized with below panels.

Input, a pair of Hi-C contact maps to be compared.

Model selection and p-value calculation.

Resulting significance map.

As can be seen on last panel our model has the abillity to enhance long range differential interactions.

5 Detection of differential significantly interacting regions

Given differential (significance) map often of interest are groups of depleted/enriched cells, which are close neighbors (in terms of euclidean distance in i, j, p-value space) for some p-value cutoff. To find such aggregates of cells we designed a simple algorithm based on connected components search. Below panel with pseudocode explains main steps of this procedure.

Selection of significantly depleted/enriched regions.

Briefly an algorithm first selects non zero cells from significance matrix, sorts them and divides into 2 categories: negative (depleted) and positive (enriched). Then for both categories a bilinear model is fitted. An intersection point between linear functions derived from the model indicate threshold (red vertical line), which further divides each category on significant and non significant differentially interacting regions (cells).

Selection of significantly depleted regions, bilinear model fitting.

Selection of significantly enriched regions, bilinear model fitting.

In the last step a zero matrix with the same shape as initial significance matrix is constructed. Cells, which coordinates correspond to that of significant ones are marked with 1. Resulting binary matrix is scanned for connected components. To perform connected components search clump function from raster package is used (Hijmans 2018). An algorithm outputs 2 elements:

list with dataframes containing connected components (rows correspond to cells within the same component),
dataframe with rectangle like regions enclosing connected components.

Providing such output we allow user flexibility to individually filter resulting regions by the size of connected components (i.e. number of significant cells within component). An example output illustrated on significance map is presented below for 2 thresholds: components with at least 2 cells and components containing at least 10 cells.

Long range differential interactions with connected components containing at least 5 significant cells.

Long range differential interactions with connected components containing at least 10 significant cells.

References

Dixon, Jesse R, Inkyung Jung, Siddarth Selvaraj, Yin Shen, Jessica E Antosiewicz-Bourget, Ah Young Lee, Zhen Ye, et al. 2015. “Chromatin Architecture Reorganization During Stem Cell Differentiation.” Nature 518 (7539). Nature Publishing Group: 331.

Hijmans, Robert J. 2018. Raster: Geographic Data Analysis and Modeling. https://CRAN.R-project.org/package=raster.

Hofert, Marius, Ivan Kojadinovic, Martin Maechler, and Jun Yan. 2017. Copula: Multivariate Dependence with Copulas. https://CRAN.R-project.org/package=copula.

Imakaev, Maxim, Geoffrey Fudenberg, Rachel Patton McCord, Natalia Naumova, Anton Goloborodko, Bryan R Lajoie, Job Dekker, and Leonid A Mirny. 2012. “Iterative Correction of Hi-c Data Reveals Hallmarks of Chromosome Organization.” Nature Methods 9 (10). Nature Publishing Group: 999.

Ivan Kojadinovic, and Jun Yan. 2010. “Modeling Multivariate Distributions with Continuous Margins Using the copula R Package.” Journal of Statistical Software 34 (9): 1–20. http://www.jstatsoft.org/v34/i09/.

Jun Yan. 2007. “Enjoy the Joy of Copulas: With a Package copula.” Journal of Statistical Software 21 (4): 1–21. http://www.jstatsoft.org/v21/i04/.

Lajoie, Bryan R, Job Dekker, and Noam Kaplan. 2015. “The Hitchhiker’s Guide to Hi-c Analysis: Practical Guidelines.” Methods 72. Elsevier: 65–75.

Le Dily, François, Davide Baù, Andy Pohl, Guillermo P Vicent, François Serra, Daniel Soronellas, Giancarlo Castellano, et al. 2014. “Distinct Structural Transitions of Chromatin Topological Domains Correlate with Coordinated Hormone-Induced Gene Regulation.” Genes & Development 28 (19). Cold Spring Harbor Lab: 2151–62.

Lieberman-Aiden, Erez, Nynke L Van Berkum, Louise Williams, Maxim Imakaev, Tobias Ragoczy, Agnes Telling, Ido Amit, et al. 2009. “Comprehensive Mapping of Long-Range Interactions Reveals Folding Principles of the Human Genome.” Science 326 (5950). American Association for the Advancement of Science: 289–93.

Lun, Aaron TL, and Gordon K Smyth. 2015. “DiffHic: A Bioconductor Package to Detect Differential Genomic Interactions in Hi-c Data.” BMC Bioinformatics 16 (1). BioMed Central: 258.

Marius Hofert, and Martin Mächler. 2011. “Nested Archimedean Copulas Meet R: The nacopula Package.” Journal of Statistical Software 39 (9): 1–20. http://www.jstatsoft.org/v39/i09/.

Schepsmeier, Ulf, Jakob Stoeber, Eike Christian Brechmann, Benedikt Graeler, Thomas Nagler, and Tobias Erhardt. 2018. VineCopula: Statistical Inference of Vine Copulas. https://CRAN.R-project.org/package=VineCopula.

Trivedi, Pravin K, David M Zimmer, and others. 2007. “Copula Modeling: An Introduction for Practitioners.” Foundations and Trends in Econometrics 1 (1). Now Publishers, Inc.: 1–111.

Yaffe, Eitan, and Amos Tanay. 2011. “Probabilistic Modeling of Hi-c Contact Maps Eliminates Systematic Biases to Characterize Global Chromosomal Architecture.” Nature Genetics 43 (11). Nature Publishing Group: 1059.

Yang, Tao, Feipeng Zhang, Galip Gurkan Yardimci, Fan Song, Ross C Hardison, William Stafford Noble, Feng Yue, and Qunhua Li. 2017. “HiCRep: Assessing the Reproducibility of Hi-c Data Using a Stratum-Adjusted Correlation Coefficient.” Genome Research. Cold Spring Harbor Lab, gr–220640.

CopulaHiC

Rafal Zaborowski

2018-12-01